Smoothed Particle Hydrodynamics in Thermal 
Phases of a One Dimensional Molecular Cloud 
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Abstract — We present an investigation on effect of the ion- 
neutral (or ambipolar) diffusion heating rate on thermal phases 
of a molecular cloud. We use the modeling of ambipolar diffusion 
with two-fluid smoothed particle hydrodynamics, as discussed by 
Nejad-Asghar & Molteni. We take into account the ambipolar 
drift heating rate on the net cooling function of the molecular 
clouds, and we investigate the thermal phases in a self-gravitating 
magnetized one dimensional slab. The results show that the 
isobaric thermal instability criterion is satisfied in the outer parts 
of the cloud, thus, these regions are thermally unstable while the 
inner part is stable. This feature may be responsible for the planet 
formation in the outer parts of a collapsing molecular cloud 
and/or may also be relevant for the formation of star forming 
dense cores in the clumps. 

I. Introduction 

In most plasmas, the forces acting on the ions are different 
from those acting on the electrons, so naively one would 
expect one species to be transported faster than the other, 
by diffusion or convection or some other process. If such 
differential transport has a divergence, then it will result in 
a change of the charge density, which will in turn create 
an electric field that will alter the transport of one or both 
species in such a way that they become equal. In astrophysics, 
ambipolar diffusion refers specifically to the decoupling of 
neutral particles from plasma. The neutral particles in this case 
are mostly hydrogen molecules in a cloud, and the plasma is 
composed of ions (mostly protons) and electrons, which are 
tied to the interstellar magnetic field. 

Nejad-Asghar (2007) has recently made the assumption 
that the molecular cloud is initially an uniform ensemble 
which then fragments due to thermal instability. He finds that 
ambipolar drift heating is inversely proportional to density and 
its value, in outer parts of the cloud, can be significantly 
larger than the average heating rates of cosmic rays and 
turbulent motions. His results show that the isobaric thermal 
instability can occur in outer regions of the cloud. The rapid 
growth of thermal instability results a strong density imbalance 
between the cloud and the low-density surroundings; therefore 
it may produce the cloud fragmentation and formation of the 
condensations. 

Many authors have developed computer codes that attempt 
to model ambipolar diffusion. Black & Scott (1982) used 
a two-dimensional, deformable-grid algorithm to follow the 
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collapse of isothermal, non-rotating magnetized cloud. The 
three-dimensional work of Mac Low et al. (1995) treats the 
two-fluid model in a version of the ZEUS magnetohydrody- 
namic code. An algorithm capable of using the smoothed 
particle hydrodynamics (SPH) to implement the ambipolar 
diffusion in a fully three-dimensional, self-gravitating system 
was developed by Hosking & Whit worth (2004). They de- 
scribed the SPH implementation of two-fluid technique that 
was tested by modeling the evolution of a dense core, which 
is initially thermally supercritical but magnetically subcritical. 
Nejad-Asghar & Molteni (2007) have recently optimized the 
two-fluid SPH implementation to test the pioneer works on 
the behavior of the ambipolar diffusion in an isothermal self- 
gravitating molecular layer. 

In this paper, we include the SPH equivalent of the energy 
equation, which follows from the first law of thermodynamics. 
Here, we include the ambipolar drift heating rate in the net 
cooling function, and we investigate the thermal phases in the 
self-gravitating magnetized molecular layer. For this purpose, 
the two-fluid SPH technique and its algorithm are given in 
section 2. Section 3 devotes to the basic equations of the one 
dimensional molecular layer and its evolution by SPH. The 
summary and conclusion are presented in section 4. 

II. Numerical method 

In two-fluid SPH technique of Hosking & Whitworth 
(2004), the initial SPH particles represent by two sets of 
molecular particles: magnetized ion SPH particles and non- 
magnetized neutral SPH particles. In this method, for each 
SPH particle we must create two separate neighbor lists: one 
for neighbors of the same species and another for those of 
different species. Consequently, each particle has two different 
smoothing lengths. In the following sections we refer to neutral 
particles as a and (3, and ion particles as a and b; the subscripts 
1 and 2 refer to both ions and neutral particles. We adopt the 
usual smoothing spline-based kernel (Monaghan & Lattanzio 
1985) and apply the symmetrized form proposed by Hernquist 
& Katz (1989) 



W 12 =W{\r 1 -r 2 \,h 12 ) 



(1) 



where r l5 r 2 are positions of the particles 1 and 2, respectively. 
h\2 is the smoothing length of particle 1 when considering 
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neighbors of kind 2. We update the adaptive smoothing length 
according to the octal-tree based nearest neighbor search 
algorithm (NNS). 

The exact fraction of the total fluid that is ionized depends 
upon many factors (e.g. the neutral density, the cosmic ray 
ionization rate, how efficiently ionized metals are depleted 
on to dust grains). Here, we use the expression employed by 
Fiedler & Mouschovias (1992), which states that for 10 s < 
n < 10 15 m~ 3 , 

Pi =e{pV* + e'p-*), (2) 

where in standard ionized equilibrium state, e ~ 7.5 x 
lO-^kg 1 / 2 .!!!- 3 / 2 and e' - 4x 10" 44 kg 5 / 2 . m - 15 / 2 are valid. 
In reality, the gas in this case is very weakly ionized, thus, we 
adopt the approximation p = p n + Pi ~ pn in fluid equations. 
In this case, the molecular cloud is considered as global 
neutral which consists of a mixture of atomic and molecular 
hydrogen (with mass fraction X), helium (with mass fraction 
Y), and traces of CO and other rare molecules, thus, the mean 
molecular weight is given by 1/p = X/2 + Y/A. 

The neutral density in place of neutral particles is estimated 
via usual summation over neighboring neutral particles 



Pn 







(3) 



while in place of ions, p n>a , is given by interpolation technique 
from the values of the nearest neighbors. The ion density is 
evaluated via equation (f2]) as follows 

p iA = 1.8 x 10- 9 p^i 2 (l + 3.5 x 10- 17 p~ 5 / 2 ), (4) 

for both places of ions and neutral particles. According to this 
new ion density, we update the mass of ion a as 

n new 



PS 



Id 



(5) 



in each time step so that the usual summation law 

Pt,a^^m b W ab , (6) 



for ions might being appropriate. 

The SPH form of the drift velocity of ion particle a is given 
by Hosking & Whitworth (2004) as 

1 r 1 

Vd,a 



"fADPn,a 



B b B a — 

, Pi,b dz a 



- mbUab 



dW ab 

dz a 



(7) 



where II a j, is the artificial viscosity between ion particles a 
and b 



n„ 



Pab 

o, 



if v ab .r afc < 0, 
otherwise, 



(8) 



where p ab — \{p a + pb) is an average density, a* ~ 1 and 
(3* ~ 2 are the coefficients, a a b = 5(^0 + a&) is the mean 
sound speed, and p a b is defined as its usual form 

Vab ' r ab 1 



Pab 



Kb rijhi b + omi 



(9) 



with h a b = h(h a + hb). Nejad-Asghar, Khesali & Soltani 
(2008) has recently considered the coefficients in the Mon- 
aghan's standard artificial viscosity as time variable, and a re- 
striction on them is proposed such that avoiding the undesired 
effects in the subsonic regions. Here, we use the Monaghan's 
standard artificial viscosity, since the cloud contraction is 
quasi-hydrostatic and there is not supersonic motions and 
shock formation during this contraction. 

It is better to keep in mind the second golden rule of SPH 
which is to rewrite formulae with the density inside operators 
(Monaghan 1992). In this case, the drift velocity is interpolated 
as 



lADPn,a 



1 y m b fz> . 



^dWab 



Pi,b 



~Pi,a > J Hab 

u Pi,b 



dz a 
dW ah 

dz a 



(10) 



where two extra density terms are introduced, one outside 
and one inside the summation sign. This comes as a result of 
the approximation to the volume integral needed to perform 
function interpolation (Nejad-Asghar & Molteni 2007). 

There is no analytical expression that allows to calculate the 
value of drift velocity of the neutral particles. Here, we use 
the interpolation technique that starts at the nearest neighbor, 
then add a sequence of decreasing corrections, as information 
from other neighbors is incorporated (e.g., Press et al. 1992). 
These drift velocities at neutral places are used to estimate the 
drag acceleration 



^drag^a / AD Pi,a^> 'd,a; 



(ID 



instead the method of Hosking & Whitworth (2004) who used 
the expression of Monaghan & Kocharayan (1995). In the 
usual symmetric form, the self-gravitating SPH acceleration 
equation for neutral particle a is 

dv a \ - /Pa Pft -t-j ,dW a p 







Pi 



(12) 



where g a is the gravitational acceleration of particle a. 

The ion momentum equation assuming instantaneous veloc 
ity update so that we have 

_ ^ m (3 







P0 



-V0W a f3 + V dta 



(13) 



where the first term on the right-hand side gives the neutral 
velocity field at the ion particle a, calculated using a standard 
SPH approximation. Finally, the magnetic induction equation 
in SPH is 

dB a _ ^ m b n dWab 
~~dT ~ 



Pb 



dz a 



(14) 



where the usual notation is used. 

In ambipolar diffusion process, the ion particles are physi- 
cally diffused through the neutral fluid, thus, the ions will be 
bared in the boundary regions of the cloud (i.e. without any 
neutral particles in their neighbors). We check the position 
of ions before making a tree and nearest neighbor search, 
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THERMALAYERSET 
set inital values 




INPUT 

use the saved values 



THERMALAYEREND 
check the time for stopping the program and/or saving data at OUTPUT 
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TREENNSL 
make tree, find neighbors and smoothing lengths 



DENSITY 

find smoothed density of neutrals in positions of neurtal and ions, and find the 
ion's density via its relation to neutral's density 



find temperature, pressure, and sound speed 



DRIFTVEL 
find drift velocity of ions 



RATES 

find time-rates of neutral's velocity, energy, and megnetic field 



update the velocity of ions via its drift velocity 



COURANT 
find minimum time-step 



advance the particles for one time-step 



X 



Fig. 1. Algorithm of the two-fluid smoothed particle hydrodynamics for 
simulation of thermal phases in a one dimensional molecular cloud. 



so that we send out the bared ion particles in the boundary 
regions of the simulation at next time-step. Both the cloud and 
boundary regions contain ion and neutral particles. We set up 
boundary particles (Ahi up and down in z) using the linear 
extrapolation approach (from the values of the inner particles) 
to attribute the appropriate drift velocity, drag acceleration, 
pressure acceleration, and the magnetic induction rate to the 
boundary particles. The algorithm of the two-fluid SPH for 
simulation of thermal phases in a one dimensional molecular 
cloud is shown in Fig. [T] 

III. LAYER EVOLUTION 
A. basic equations 

We write the continuity equation of neutral part as its 



common form 



dp 



dv 



dt p dz 



(15) 



while, the relation (O is used to determine the ion density 
wherever it is required in the two-fluid equations. The mo- 
mentum equation of neutral part then becomes 

dv 19, B 2 

17=9--^r(P+T^—) (16) 
at p oz 2po 

where the gravitational acceleration g obeys the poisson's 
equation 

(17) 



oz 



and the pressure is given by the ideal gas equation of state 



P = 



pm H 



-pT 



(18) 



where TOh is the mass of hydrogen atom and the temperature 
T is approximated the same as for both neutral and ion fluids 
(T t = T n = T). 

The thermal energy per unit mass is 



(19) 



where the mean internal (rotation and vibration) energy of an 
H2 molecule is included. The energy equation follows from 
the first law of thermodynamics, that is 

du 
~dl 



pdv 

p az ' 



(20) 



where ^( Pi t) is the net cooling function 

Q (p,t) = A w(t ^)^ ( " ) " + r AD), (21) 

where Tcr and Tad are the heating rates due to cosmic rays 
and ambipolar diffusion, respectively, and A(„j and /3( ra ) are 
the parameters for the gas cooling function that here we use the 
polynomial fitting functions, outlined by Nejad-Asghar (2007), 
as follows 

log ( 7T~T -1 ) = -8-98 - 0.87(log ^) 



3.07-0.11(log— ) 

n 



n 



-0.14(log- 



n 

0.13(log-) 2 

n 



(22) 
(23) 



where no — 10 m -3 . 

The magnetic fields are directly evolved by charged fluid 
component, as follows: 

d -(Bv d ), (24) 



dB _ B dv 
dt dz 



dz 



where the last term outlines the ambipolar diffusion effect with 
drift velocity, 



1 



r) R 2 

v d = £(5-), (25) 

jADeppi oz 2p,o 

which is obtained by assumption that the pressure and grav- 
itational force on the charged fluid component are negligible 
compared to the Lorentz force because of the low ionization 
fraction. The notation jad ~ 3.5 x 10 10 m 3 /kg.s in drift 
velocity represents the collision drag coefficient. 

B. results 

The chosen physical scales for length and time are [/] = 
200AU, and [t] = 10 3 yr, respectively, so that velocity unit 
is approximately [v] = lkm.s -1 . The Newtonian constant of 
gravitation is set G = 1[to]~ 1 [Z] 3 [^] -2 for which the calculated 
mass unit is [to] = 4.5 x 10 29 kg. Consequently, the derived 
physical scale for density, energy per unit mass, and drag co- 
efficient are [p] = 1.7 x 10~ n kg.rrr 3 , [u] = 10 6 J.kg~\ and 
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Iad = 1-8 x 10 10 |7] 3 [m] _1 respectively. In this manner, 
the numerical values of e and e' are 1.8 x lO^ 9 ^] -3 / 2 ^] 1 / 2 
and 3.5xl0~ 17 [Z]~ 15 / 2 [m] 5 / 2 , respectively. The magnetic field 
is scaled in units such that the constant fio is unity. Since the 
magnetic flux density has dimensions 



while /i has dimensions 



[t] [charge] 

MM 



(26) 



(27) 



[charge] 2 ' 

specifying po — 1 therefore scales the magnetic field equal to 
[B] = 5.1nT. With aforementioned units, the thermal energy 
per unit mass (19[ is represented by 8.3 x 10 _3 (5A"/4 + 
3Y/8)T, the heating rates due to cosmic rays and ambipolar 
diffusion are T C r = 1.1 x 10~ 3 [u] /[i] and 

^AD,a — "fADPi,aV^ a , (28) 

respectively, and the parameters for the gas cooling function 



l"MT^)=-4-48-0.87(log 
-0.14(log 



2.24 x 10 



Pa 



2.24 x 10 



3l) 2 , (29) 



P a = 3.07-0.11(log 
-0.13(log 



Pa 



2.24 x 10- 4 

Pa \2 



-r 



2.24 x 10-- (30) 
The initial conditions for this simulation are a parallel 
magnetic field directed perpendicular to the z-axis so that 
the initial ratio of magnetic to gas pressure is everywhere a 
constant (an. = 1)> an d a density profile given by 

P = 2 (31) 

cosh (Z/Zoo) 

where p is the initial central density and z^ = 
c s \J(\ + ao)/27rG j oo is a length-scale parameter according 
to the initial sound speed c s (see Fig. 1 of Nejad-Asghar & 
Molteni 2007). The magnetic field is assumed to be frozen in 
the fluid of charged particles and the central density is assumed 
to be po = 2.24 x 10~ 4 [p]. We choose a molecular cloud 
with the mass fraction of molecular hydrogen and helium are 
X = 0.75 and Y = 0.25, respectively, and have initial uniform 
temperature of To = 50K. We assume that the cloud layer 
is spread from z — — 93.2[Z] to z — +93.2 [Z] (according to 
Fig. 1 of Nejad-Asghar & Molteni 2007). The initial values 
of the cooling and heating functions are shown in Figure |2^. 
As presented in this figure, the isobaric thermal instability 
criterion, 

op op 

is satisfied in the outer parts of the cloud, thus, these regions 
are thermally unstable while the inner part is stable as outlined 
by Nejad-Asghar (2007). 



The present SPH code has the main features of the TreeSPH 
class so that the nearest neighbors searching are calculated 
by means of this procedure. We integrate the SPH equations 
using the one order simple integration scheme. The selection 
of time-step, At, is of great importance. There are several 
time-scales that can be defined locally in the system. For each 
particle 1, we calculate the smallest of these time-scales using 
its smallest smoothing length, h\, i.e. 

. . hi hi hi 



AU 



C a 



\vi\ ' V A ,l ' Ol 



(33) 



where va — B/^/p pi is the Alfven speed of ion fluid and 
C cour is the Courant number which in this paper is adopted 
equal to 0.3 (for numerical stability). The evolution were 
carried out to time < u > / < A, Y >~ 0.43/0.025 = 17.2[i] 
so that the stability of the outer parts of the cloud may be 
revealed. The values of the cooling and heating rates versus 
position (and neutral density) are shown in Fig. |2j5-d, at times 
* = 3.5 [t], 10.5 [t] and 17.2 [t], respectively. 

IV. Summary and conclusions 

In the preceding work of Nejad-Asghar & Molteni (2007), 
two-fluid SPH implementation was used to check the pioneer 
works on behavior of the ambipolar diffusion in an isothermal 
self-gravitating layer. In this paper, we include the SPH equiv- 
alent of the energy equation that can be obtained from the first 
law of thermodynamics. Here, we incorporate the ambipolar 
drift heating in the net cooling function, and we investigate the 
thermal phases in the self-gravitating magnetized molecular 
layer. 

The values of the cooling and heating functions are shown 
in Figure |2] According to this figure, the isobaric thermal 
instability criterion is satisfied in the outer parts of the cloud, 
thus, these regions are thermally unstable while the inner part 
is stable. The SPH equations were integrated, and the evolution 
were carried out to time 17.2[i] so that the stability of the outer 
parts of the cloud is revealed. 

It is obvious that the instability of the cloud at the outer 
parts, causes the formation two relative cool regions in that 
area. The rapid growth of thermal instability results a strong 
density imbalance between the cloud and the low-density 
surroundings. This feature may be responsible for the planet 
formation in the outer parts of a collapsing molecular cloud 
and/or may also be conscientious for the formation of star 
forming dense cores in the clumps. 
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Fig. 2. The values of the cooling (solid) and heating (dash) functions versus position (and neutral density), at initial time and t = 3.5[t], 10.5[t] and 17.2[t], 
respectively. Thermal unstable regions are shaded in this figure. 
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